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Abstract 

We consider two qubits interacting with local and collective thermal reservoirs. 
Each spin-reservoir interaction consists of an energy exchange and an energy con- 
serving channel. We prove a resonance representation of the reduced dynamics 
of the spins, valid for all times t > 0, with errors (small interaction) estimated 
rigorously, uniformly in time. Subspaces associated to non-interacting energy dif- 
ferences evolve independently, partitioning the reduced density matrix into dynam- 
ically decoupled clusters of jointly evolving matrix elements. Within each subspace 
the dynamics is markovian with a generator determined entirely by the resonance 
data of the full Hamiltonian. Based on the resonance representation we examine 
the evolution of entanglement (concurrence). We show that, whenever thermaliza- 
tion takes place, entanglement of any initial state dies out in a finite time and will 
not return. For a concrete class of initially entangled spin states we find explicit 
bounds on entanglement survival and death times in terms of the initial state and 
the resonance data. 



1 Introduction and outline of main results 

We consider two qubits Si, S2 (two spins 5) in a thermal environment. If the spatial 
separation of the qubits is small compared to the correlation length of the environ- 
ment then the qubits feel the same interaction with the reservoir, called a collective 
interaction. For separated qubits the interaction is described by independent reservoirs 
coupled to each qubit, called a local interaction. In experiments both kinds of inter- 
action occur at the same time and so we consider three independent heat reservoirs, 
Ri , R2 (local) and R (collective) . 

Each of the interactions Si + Ri, S2 + R2 and (Si, S2) + R has two channels given 
by energy conserving and energy exchange couplings. The former are also called "non- 
demolition" interactions, leaving the energy of Si and S2 invariant. Under this evolution 
the populations are constant (diagonal elements of reduced density matrix of Si + S2 
in the energy basis). Nevertheless, correlations decay in time (off-diagonal density 
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matrix elements). This process is called "phase decoherence" . Models with constant 
populations have a multitude of invariant states (at least all diagonal density matrices) . 
They cannot describe asymptotic processes such as the approach to a final (equilibrium) 
state. Those are induced by energy exchange interactions. This is why both interaction 
channels should be considered simultaneously. 

Our main results are Theorems 2.1 and 2.3 on the reduced dynamics, and Theorem 
2.4 on entanglement evolution. In the former we prove a "resonance approximaton" of 
the reduced dynamics with a remainder small in the interaction, valid uniformly in time 
t > 0. Theorem 2.1 shows that the dynamics is expressed by just a few parameters (the 
resonance data). In Theorem 2.3 (and Section 4.1) we calculate explicitly the resonance 
data for the model of two qubits interacting with local and collective reservoirs. We 
use these results to establish, in Theorem 2.4, rigorous bounds on the times of survival 
and death of entanglement for a class of initial qubit states. 

Reduced dynamics. The dynamics of the qubits is obtained by tracing out the 
reservoirs degrees of freedom, and is given by a time-dependent reduced density matrix 
pt. Assuming that the qubits are initially not entangled with the reservoirs (but may 
be entangled among themselves), the initial qubit state po evolves according to a linear 
mapping of density matrices [10] 

Po ^ Pt = V x (t)p . 

In this introduction we consider the coupling strenghts of all interactions to be propor- 
tional to some overall coupling constant x. 

The dynamical map V x (t) is not a (semi-) group in t, but it is common in physics 
to make a so-called markovian master equation approximation (Born-, Markov- and 
rotating wave approximations), which restores the group property: V x {t) « e tc>! . Here, 
C x = Co + x 2 K^ is a Lindblad generator which includes the uncoupled dynamics and 
a second-order correction term (traditionally denoted by if" (Davies)). The validity 
of this approximation is based on physical considerations and is argued to work well 
if the bath correlation times are much smaller than the system relaxation time [10]. 
However, to our knowledge, no rigorous control of the error has been achieved in this 
method. In the extensive physics literature one finds different forms of (derived 
heuristically), only one of them - the so-called Davies generator [11, 12] - giving a 
positivity-preserving dynamics [14] (guaranteeing that e tC "po is a density matrix for 
all t > 0). The mathematical procedure leading to this correct generator is the so-called 
weak coupling-, or Van Hove limit. It is stated as follows [11, 12, 19, 13]: for any a > 0, 

lim sup \\V x (t) - e tC "\\ = 0. 

(It is assumed that p acts on a finite dimensional Hilbert space.) The disadvantage 
of this approach is that one can only consider time scales up to 0{x~ 2 ). In [15, 16] 
the weak coupling analysis is extended to time scales up to 0(x -4 ) by replacing the 
semigroup e tl ~ by a more complicated mapping taking into account non-Markovian 
effects. Still, in order to examine large times, one has to diminish simultaneously the 
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value of the coupling constant. We prove in Theorem 2.1 a resonance approximation, 
improving the weak-coupling result to 

sup || V*(i) -e* M "|| < C7x 2 , 

where M H = Cq + >c 2 K^ + ... is analytic in x. The gain is control uniform in t > 0, it 
is achieved by replacing C x by the more complicated generator M K which contains all 
orders of x. M K commutes with the free generator Cq, and the second-order correction 
is the Davies generator. The latter is linked to quantum resonance theory by the 
relation (Proposition 7.1) 

K*= (iA e )*. 
eespec(£o) 

Here, (A e )* is the adjoint operator of the level shift operator A e (acting on ^(T-Ls) = 
Us <8> T~Ls, see (3-5)), whose spectrum yields the complex energy corrections (second 
order) of the unperturbed Bohr energy e£l (energy difference of the system Hamil- 
tonian). The eigenvalues of M K are the resonances (complex eigenvalues) of a non- 
selfadjoint Liouville operator (see (3.1)). 

In the master equation setting, the rotating wave (or secular) approximation con- 
sists in neglecting, in the evolution equation, quickly oscillating terms proportional to 
e i(e-e )t if e ^ where e, e' are Bohr energies [10]. This leads to an approximate 
dynamics in which spectral subspaces associated to different Bohr energies evolve inde- 
pendently. The fact that M x leaves eigenspaces of Cq invariant gives thus a proof of the 
validity of the rotating wave approximation. This decoupling of evolution of spectral 
subspaces may be very useful in the analysis of open systems with many degrees of 
freedom, as the density matrix is partitioned into independent clusters of jointly evolv- 
ing matrix elements. Some clusters may decay quickly and one is soon left with only a 
'sparse' density matrix. Also, for certain applications (quantum protocols) only a few 
clusters may be important, and one can focus on their analysis directly by studying 
the corresponding level shift operators. 

The present resonance approach builds on [23, 24, 25], giving the dynamics of 
reduced density matrix elements directly in terms of spectral data. An improvement 
of the weak coupling limit to all times for a single spin coupled to a Bosonic reservoir 
has also given in [19]. 

Evolution of entanglement. Entanglement is a central notion of modern quan- 
tum theory and particularly quantum information and computation. It is a measure 
for quantum correlations between subsystems. In the simplest setting, entanglement 
measures how far away from being a product state a given state of a bipartite sys- 
tem is. (A product state also being called disentangled.) There are various notions of 
entanglement [18]. We concentrate on entanglement of formation [7], defined for the 
density matrix p of a bipartite system A + B to be 

£{p)= inf V^-5(Tr B |^)(^|), 

where S(x) = — Tr(xlnx) > is the von Neumann entropy (when dealing with systems 
of spins \ it is natural to take the binary logarithm). Here, Tr# is the partial trace 



3 



over system B, and the infimum is taken over all probabilities < pj < 1, YljPj = 1j 
and all vector states tpj of A + B, \\tpj\\ = 1, such that p = X^PjlV'jXV'jl- 

The expression for entanglement of formation involves a hard problem of minimiza- 
tion over all possible realizations of a given state p, and its calculation is a (too) difficult 
enterprise. However, Wootters [29] has shown that if both A and B are spins ^ (each 
having state space C 2 ), then £{p) is related strictly monotonically to the concurrence 
C(p) defined by an expression obtainable directly from the density matrix p, see (2.30)- 
(2.33). The concurrence of two spins satisfies < C(p) < 1. It vanishes if and only if 
p can be written as a mixture (convex combination) of pure product states. 

The link between entanglement and concurrence established by Wootters has pro- 
voked a wealth of investigations on entanglement of two spins in interaction with noises 
(classical or quantum). In particular, the phenomena of entanglement decay, sudden 
entanglement death, death and revival as well as creation of entanglement due to a col- 
lective reservoir have been discovered. There is a large body of work also concerning 
entanglement of continuous variable systems (bosonic modes), using different entan- 
glement measures, such as "negativities". We refer to the review [1] and the extensive 
bibliography therein. In the present work, we focus on two qubits and concurrence as 
a measure of entanglement. 

Other than numerically, entanglement has been studied in the weak coupling marko- 
vian master equation regime (Lindblad dynamics) or for explicitly solvable models. 
Within the context of the markovian description, once the Lindblad dynamics is as- 
sumed, a mathematically rigorous analysis is sometimes be possible. Our goal here is 
to start with the full microscopic (Hamiltonian) description and extract the dynamics 
of entanglement rigorously, without assuming a Lindblad evolution as a starting point. 
We next point out some references on entanglement evolution; the literature on the 
topic is huge and many more references are found in the papers cited below. 

In the markovian approximation, Yu and Eberly [30] find an upper bound on decay 
of concurrence of two qubits: C(pt) < e _7 *C(po) for some 7 > 0. (Each spin is coupled 
locally to a zero-temperature cavity reservoir through an energy exchange interaction.) 
They show that some initial states po satisfy C(pt) = for all times exceeding an 
entanglement death time, but also that there are initial states for which C{pt) > 
for all times. Furthermore, in some models with purely energy-conserving interaction, 
they find entanglement-free subspaces [31]. (Explicitly solvable model with classical 
(commutative, stochastic) local and collective noises.) 

Bellomo et al. [4] consider two spins locally interacting with zero temperature cav- 
ities in a non-markovian regime (explictly solvable energy-exchange interaction). They 
show that entanglement of certain initial states undergoes sudden death and revival: 
initial entanglement decays to zero, may stay zero for a while, but then reappears (with 
some loss), and so on. The interpretation is that initial entanglement is shifted from 
the spins to the reservoirs (intially unentangled), then shifted back to the spins with 
some loss. 

Braun [9] considers two spins coupled collectively to a single thermal reservoir 
(energy-conserving interaction, explicitly solvable). He shows that certain initially 
unentangled states of the spins will acquire some entanglement for a while, due to the 
interaction with the common reservoir. In a similar spirit, using the Peres-Horodecki 
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(partial transposition) criterion for entanglement detection [27, 18], Benatti et al. show 
that entanglement can be created by a collective environment for certain initial condi- 
tions in markovian, not explicitly solvable models [5, 6]. 

Our present contribution are rigorous bounds on entanglement survival and death 
times. We show in Theorem 2.2 that the entanglement of any initial state will decay 
in a finite time if the entire system has the property of return to equilibrium. This 
property holds generically, and so the question is: how long can initial concurrence 
prevail, and when will it certainly have died out without returning? In Theorem 2.4 we 
give bounds on entanglement survival and death times, linking them to the resonance 
data for the model with all interactions present. To our knowledge, this is the first 
time rigorously established bounds are given for not explicitly solvable models. 

We have announced Theorems 2.1, 2.3 and 2.4 in [21] without proofs. The main 
result of that paper is a numerical analysis showing that the resonance approximation 
captures the phenomena of sudden death and revival, and of creation of entanglement. 

2 Main results 
2.1 Model 

The space of pure states of S = Si + S2 is given by 

^ S = C 2 ®C 2 , (2.1) 

its Hamiltonian is 

H S = BrfZ + BzSZ. (2.2) 

Here B\,B2 > are the values of an (effective) magnetic field at the positions of the 
two spins, 

" 1 



s z 



-1 



is a Pauli matrix and Sf = S z ® 1 C 2, 5| = He 2 ® S z . 

Each thermal reservoir is described by a spatially infintely extended free bose gas in 
equilibrium at temperature 1//3 > 0, 1 having positive particle density (phase without 
Bose-Einstein condensate). A mathematical description of reservoirs of free particles 
is usually given in terms of a state of the CCR Weyl algebra [8, 2], but we explain 
our model here in terms of a*(k), a(k), the bosonic creation and annihilation operators 
of a particle with momentum fc £ i 3 , satisfying the canonical commutation relation 
[a*(k), a(l)] = S(k—l). The operators a(k), a*(k) act on bosonic Fock space (represented 
here in Fourier, or momentum space) 



oo 



J=$L^(M 3n ,d^), (2.3) 



n=0 



1 Our approach and our results can easily be applied to the situation where each of the three reservoirs 
has a different temperature, but the scope of this work is to compare local versus collective phenomena, 
and thus we take all reservoir temperatures to be the same. 



5 



where L^ yra is the space of all functions invariant under permutation of the arguments 
ki,...,k n (symmetric functions) . The Hamiltonian of the reservoir is given by 

H R = [ \k\a*(k)a(k)d 3 k, (2.4) 

where \k\ is the energy of a single particle with momentum k. The equilibrium state of 
a reservoir at temperature > is the quasi- free state uip determined by 

W/J (o(fc)) = «„(*•(*)) = 0, w/j(a*(fc)a(Z)) = (2.5) 

The second relation implies that the probability density distribution of particles with 
momentum k is given by (e^' fc ' — l) -1 , which is with Planck's black-body radiation law. 
For g € L 2 (M 3 , d 3 fc) we define the smoothed-out creation and annihilation operators 

a* (g) = [ </(£>* (fc)d 3 fc, a(g) = f gjk)a(k)d 3 k, 



where g(k) is the complex conjugate of g{k) and we define the field operator 

<t>(9) = ^[a*(g)+a(g)}. 
The Hilbert space of the total system is given by the tensor product 

n = n s ® Jki ® ^r 2 ® Jk, (2.6) 

where each F Rl , J 7 ^ , F R is a copy of J 7 , and %s is given in (2.1). The full Hamiltonian 
takes the form 

H = H s + H Rl +H R2 + H R (2.7) 

+ (AiSf + A 2 5f ) ® 0( 5 ) (2.8) 

+ ( Kl 5f + «25|) ® <K/) (2.9) 

+ /xi5f ®^i) + /x 2 5f (8)0( 5 2 ) (2.10) 

+ ^Sf®0(/i) + i*Sf ®#/ 2 ). (2.11) 

The first four terms, (2.7), are the uncoupled Hamiltonians, where H$ is given in (2.2) 
and each of the H Rl , H R2 , H R is the operator (2.4), acting on the appropriate factor of 
the Hilbert space (2.6). 

The collective interaction is determined by the operators (2.8) and (2.9), where 4>{g) 
acts nontrivially only on the factor .Fr in (2.6). The terms (2.8) describe collective 
energy exchange interactions of the spins with R, with respective coupling constants 
Ai, A2 G M. Energy exchanges are implemented by the spin flip operator (Pauli matrix) 

b ~ 1 ' 
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acting on either of the factors C 2 of Hs (denoted by Sf or S% , analogously to Sf, Sf in 
(2.2)). The operator (2.9) commutes with H$ and describes collective energy preserving 
interactions, with respective coupling constants k±,K2 and held operator 4>{f) acting 
nontrivially on the factor of % . 

The local interactions are governed by the operators (2.10) (local, energy exchange) 
and (2.11) (local, energy conserving). Again, fj,±, (12, v\, ^2 G ^ are coupling constants. 
The field operators 4>(gj) and 4>(fj), j = 1, 2, act nontrivially on the factor Jr. of H. 

For any state on S = Si + S2 (positive linear functional on the algebra of 
bounded operators B(Hs)) and any system observable A G B(Hs), we denote by u)g(A) 
the reduced dynamics at time t, given by 

<4(A) = uj s ® w^Ri ® c^,r 2 ® u^r (e iiH (^ ® l Rl 1r 2 ® IrK^) , (2.12) 

where is the full Hamiltonian (2.7)-(2.11). This expression is formal, and we under- 
stand that the termodynamic limit of all reservoirs is taken. 2 The system dynamics 
(2.12) determines the reduced density matrix p t by tOg(A) = Trs(ptA). 

2.2 Results 

Our main results on decoherence and disentanglement are based on a careful analysis 
of the reduced dynamics of S = Si + S2, based on a refinement of the dynamical 
theory of quantum resonances developed in [23, 24]. This theory uses analytic spectral 
deformation methods and requires the following regularity assumption. 

(A) Write h(r, S) for a function h € L 2 (R 3 ,d 3 A;) in spherical coordinates (r, S) G 
R_i_ x S 2 and consider the function 



x S> ^ { u, E) * M «, S> := |»| 1/2 { 



/i(u, £), if -u > 

if u < 



(2.13) 

where a £ I is any fixed phase, and (3 > is the inverse temperature. It 
is assumed that for h = g, g±, 32, /, fi, ji (see (2.8)-(2.11)), the function 6 >->■ 
/ig(ii + 9, S) has an analytic extension in the variable 0, as a map C x S 2 1— > 
L 2 (M x S" 2 ,dudS), from G R to £ {z £ C : < 9z < O }, such that the 
extension is continuous as 96* | 0. Here, #0 > is an arbitrary constant bounded 
above by 27r//3 (singularity of the square root in (2.13)). 

A family of form factors h satisfying this condition is h(r, S) = r p e _rm /ii(S), with 
p = — 1/2 + n, n = 0, 1, . . . and m = 1,2, and where hi is an arbitrary (integrable) 
function of the angle S. This family contains the usual physical form factors, [28]. 



2 A rigorous discussion of this point, for a general finite-dimensional system S is coupled linearly in 
the field operators is given in [23, 24]. 
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2.2.1 Dominant reduced dynamics 

Denote the eigenvalues of H$ by 

Ei = B\ + B2, E2 = B\ — B2, E3 = —B\ + B2, E4 = —B\ — B2, (2-14) 
with corresponding ordered basis {#1, . . . , #4} of Hs- Explicitly, 

#i = | + +), #2 = | + -), ^3 = |- +), #4 = 1- -), (2.15) 
where \a1a2) = \o~i) ® \<J2) and S z \±) = ±|±). Define 

x := max{|/« j |, |Aj|, |/Xj|, 1^1 : J = 1,2}. (2.16) 

Under the non-interacting dynamics (x = 0), the evolution of the reduced density 
matrix elements of S (expressed in the energy basis (2.15)) is given by 

[pt]kl = (<P k ,e- itHs po e itHs ^> = e ite <* [ Po ] k i, (2.17) 

where = E\ — E^. As the interactions with the reservoirs are turned on (some of 
Kj , Xj , fij , Vj nonzero), the dynamics (2.17) undergoes two qualitative changes. 

- The "Bohr frequencies" 

e 6 {E — E' :E,E'e spec(F s )} (2.18) 

in the exponent of (2.17) become complex resonance energies, e 1— > e e with non- 
negative imaginary parts. If Qe e > then the corresponding dynamical process 
is irreversible (decay). 

- Matrix elements do not evolve independently since the effective energy of S is 
changed due to the interaction with the reservoirs. 

Both effects are small if the coupling is small. In particular, e e — > e as x — > 0. In order 
to set up a perturbation theory we view the energy differences (2.18) as the eigenvalues 
of the Liouville operator 

L s = Hs®ts-t s ®H s , (2.19) 

acting on the doubled Hilbert space Us ® %s- Denote the spectral projection of Ls 
associated to e by P e . We have dimP e = mult(e) (multiplicity of eigenvalue). Generally, 
as the coupling parameters are turned on, there is a multitude of distinct resonance 
energies bifurcating out of e, with s = 1, . . . ,^(e) and v(e) < mult(e). The shift 
of eigenvalues e under perturbation is of order at least two in the coupling constants 
in our model. 3 The lowest order corrections 5^ are 0(x 2 ). They are the eigenvalues 
of an explicit level shift operator A e (see (3.5)), acting on RanP e : there are two bases 
{rji^} and {rji^} of RanP e , satisfying 

A e # = 6^4 S \ [A e ]*# = = S S:S , (2.20) 

3 This follows from the fact that the coupling operators (2.8)-(2.11) are linear in the field operators, 
thus having zero average in the thermal state, which implies that the first order corrections to the 
energies vanish as well. 
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Moreover, we have 

e W = e + # + 0(x 4 ) (2.21) 

and G# > 0. 

(s) 

If all the resonance energies e e become distinct at this lowest order of perturbation, 
for generic small values of the perturbation parameters, then we say that the "Fermi 
Golden Rule Condition" is satisfied. By "generic small" values of the parameters 
Kj, Xj, fJ-j,Vj, j = 1,2, we mean that they belong to a small ball in 1R 8 around the 
origin, except possibly to a set of measure zero inside the ball. (Note that the origin 
is not a generic point.) This assumption, expressed as follows, is satisfied in many 
applications (as in the model of the present paper). 

(F) For generic small values of the coupling constants, there is complete splitting of 

(s) 

resonances under perturbation at second order, i.e., all the 5e are distinct for 
fixed e and varying s. 

This condition implies that v(e) = mult(e) for all e. To quantify the clustering of 
density matrix elements into groups who evolve jointly under the full evolution, we 
introduce for each energy difference e, (2.18), the cluster set 

C(e) = {(k,l) : E k -Et = e}. (2.22) 

Denote by [pt]mn the element m, n of the reduced density matrix at time t in the 
energy basis {^j}, i.e. 

[pt]mn=0j\\^ n )^m\)- (2.23) 

The following result is obtained from a detailed analysis of a representation of the 
reduced dynamics given in [23, 24, 25]. We present this result for the specific system 
at hand, but mention that the proof (Section 3) relies entirely on generic properties 
of resonance vectors and energies, and the result can be proven for general TV-level 
systems coupled to heat reservoirs. 

Theorem 2.1 (Dominant reduced dynamics) Suppose that Conditions (A) and 
(F) hold. There is a constant xq > such that if h < kq, then we have for all 
t > 

[pt]mn= Yl A t {m,n;k,l) [p ] M + O{x 2 ), (2.24) 

(k,l)eC(E m -E n ) 

where the remainder term is uniform in t > 0, and where the amplitudes A t satisfy the 
Chapman-Kolmogorov equation 

A t+r (m,n;k,l) = ^ A t (m,n;p,q)A r (p,q;k,l), (2.25) 
(p,q)eC(E m -E n ) 

for t, r > 0, with initial condition 

A (m, n; k, I) = 5 m=k 5 n= i (2.26) 
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(Kronecker delta). Moreover, the amplitudes are given in terms of the resonance vectors 
and resonance energies by 

mult(E n -E m ) 

A t (m, n; k,l) = J2 e it£ ^-^ (<f>; ® ^VeI-eJ) {^eI-e^u ® $ m ) • 

3 = 1 

(2.27) 

Discussion. (1) The result shows that to lowest order in x, and uniformly in time, the 
reduced density matrix elements evolve in clusters. A cluster is determined by indices 
in a fixed C(e). Within each cluster the dynamics has the structure of a classical (com- 
mutative) Markov process. In general the 'transition probabilities' At(m,n;k,l) are 
complex valued. The typical stochasticity property Y2(k i)eC(E„-E m ) At( m , n \ k,l) = 1 is 
not satisfied. (Generally, due to decoherence, all matrix elements corresponding to the 
(non-diagonal) clusters decay exponentially quickly to zero at the rate min s 3f£^_^ •) 
However one sees readily that if the diagonal elements form a cluster by themselves 
(which is the case e.g. if the kernel of Ls has dimension dim'Hs)) then this cluster 
satisfies the stochasticity condition: Yl(k fe)eC(o) At(m,m; k, k) = 1. 

The clustering implies the dynamical decoupling of eigenspaces associated to differ- 
ent Bohr energies. This, together with the Chapman-Kolmogorov relation, shows that 
each subspace has a (Markov) generator of dynamics. Hence the existence of M x and 
its commutativity with Cq, as mentioned in the introduction. 

(2) The resonance energies contain, in general, terms of all (even) orders in x. 
It may happen that, due to degeneracies of energy levels of Hs, condition (F) is not 
satisfied, and that zero is still a degenerate eigenvalue at order k 2 [22]. In this situation 
a result similar to Theorem 2.1 can be proven. 

(3) The existence of an equilibrium state of the coupled system implies that one 
of the resonances e^ 1 is always zero [22], we set = 0. The condition Qte^ > for 
all e, s except e = 0, s = 1 is equivalent to the system converging to its equilibrium 
state, by which we mean that lim^oo ^{A) = ujr(A) for all observables A of S, where 
ujp is the state of S obtained by reducing the coupled equilibrium state of S plus the 
reservoirs, up coincides with the Gibbs state of S at temperature 1//3, up to 0(x 2 ). 

(4) For generic couplings, S undergoes decoherence in the energy basis, which means 
that the off-diagonal density matrix elements converge to zero for large times, at lowest 
order in the coupling. However, the diagonal elements cannot experience the same fate, 
since they must sum up to one at all times (Trpj = 1). 

- The thermalization rate is defined by 

7 th = nun %e ( s) > 0. (2.28) 

(Recall that = 0, see remark (3) above.) The cluster decoherence rate associ- 
ated to C(e), e / 0, is defined to be 

7 f ec = minSeW > (1 < s < mult(e)). (2.29) 

s 

If 7 is any of the above rates, then r = I/7 is the corresponding (thermalization, 
decoherence) time. We use the convention r = 00 <^ 7 = 0. 
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- We say the system has the property of return to equilibrium if any normal initial 
state converges to the (joint) equilibrium state as time tends to infinity. 4 



2.2.2 Finite disentanglement time 

Let p be the density matrix of two spins \. The concurrence [29] is defined by 

C(p) = max{0, D(p)}, (2.30) 

with 

and where v\ > v<i > v% > 1/4 > are the eigenvalues of the matrix 

£(p) = p(S y ® S y )p(S y ® S y ). (2.32) 

Here, p is obtained from p by representing the latter in the energy basis and then taking 
the elementwise complex conjugate, and S y is the Pauli matrix 



S y 



-i 

1 



(2.33) 



It is rather easy to see that all eigenvalues of £(p) are non-negative (but £(p) is not 
hermitian). The concurrence (2.30) takes values in the interval [0, 1]. If C(p) = then 
the state p is separable, meaning that p can be written as a mixture of pure product 
states. If C(p) = 1 then p is called maximally entangled. 

Let po be an initial state of S. The smallest to > s.t. C(pt) = for all t > to is 
called the disentanglement time (also 'entanglement sudden death time', [31, 30]). If 
C(pt) > for all t > then we set to = 00. The disentanglement time depends on the 
initial state. 

As mentioned in Section 1, it is well known that under energy-conserving interac- 
tions, initial entanglement may decay gradually to zero without being zero at any finite 
time, or it may even stay constant. However, systems with energy exchange (where 
H$ is not conserved) generically have the property of return to equilibrium [20, 3, 17]. 
A consequence of thermalization is death of entanglement (in finite time). The mech- 
anism is very simple: the equilibrium (Gibbs) state is the centre of a neighbourhood 
of disentangled states, of size oc [Tr e - ^ s ] -1 . By the property of return to equilib- 
rium, the qubit state approaches its final state which is the qubit Gibbs state plus 
a correction of the order x 2 (reduction of joint equilibrium of qubits and reservoirs). 
Thus under the condition k 1 <K [Tre _/3i ^ s ] _1 the qubit state will enter and stay in the 
entanglement free neighbourhood. We give some estimates in Section 5. 



More precisely, let j3 be the inverse temperature of the reservoirs Ri,R2,R, and consider the 
reference state Lu TC i = <^>Si ® ^s 2 ® ^Ri,/3 ®R 2 ,/3 ®^>R,/9, where ujs j are arbitrary states on Sj, j = 1, 2, 
and Wn.p is the /3-KMS state of reservoir ~Rj. Let be the algebra of observables of the total system, 
represented as a von Neumann algebra acting on the GNS representation Hilbert space H TC f of cj ro f. 
Let a* be the Heisenberg dynamics of the total system, a* is a group of *automorphisms of 9Jt. The 
system has the property of return to equilibrium if for any state u represented by a density matrix on 
Hrd and any A € 9JI, we have limt^oo w(a*(j4)) = w^ iX (j4), where up,* is the /3-KMS state w.r.t. the 
coupled dynamics of the whole system. 
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Theorem 2.2 (Finite disentanglement time) Suppose the system has the property 
of return to equilibrium at some temperature T = 1/(3 > and for some Let 
Po be any initial state of the qubits. Then there is a constant c (independent of (3 
varying in compact sets and of po) s.t. if x 2 < c [Tr e - ^ Hs ] -1 , then we have C(pt) = 
(concurrence) for all t > to, where to < oo depends on (3, x and po- 

If in addition return to equilibrium happens exponentially quickly at the rate l/r t h, 
then there is a constant d > s.t. t < r t h ln[c'Tr e _/3jHs ]. 

Theorem 2.2 shows that if the coupling constants are small, for fixed temperature, then 
the disentanglement time is finite. However, it is shown in [26] that if, at fixed coupling 
constants, the temperature is decreased sufficiently, then entanglement can persist for 
all times even if the system has the property of return to equilibrium. 



2.2.3 Decoherence and thermalization rates 

We consider the Hamiltonian Hs, (2.2), with parameters < B\ < B 2 s.t. B 2 /B\ / 2. 
These assumptions represent a non-degeneracy condition which is not essential for the 
applicability of our method, but lightens the exposition. The eigenvalues of Hs are 
given by (2.14) and the spectrum of Ls is {ei, ±e 2 , ±e3, ±e4, ±es}, with non-negative 
eigenvalues 

ei = 0, e 2 = 2Bi, e 3 = 2B 2 , e 4 = 2(B 2 - B x ), e 5 = 2{B 1 + B 2 ), (2.34) 

having multiplicities m\ = 4, m 2 = = 2, TO4 = 771,5 = 1> respectively. According to 
(2.34), the grouping of jointly evolving elements of the density matrix above and on 
the diagonal is given by five clusters 5 

d := C(ei) = {(l,l),(2,2),(3,3),(4,4)} (2.35) 

C 2 := C(e 2 ) = {(l,3),(2,4)} (2.36) 

C 3 := C(e 3 ) = {(l,2),(3,4)} (2.37) 

C 4 := C(-e 4 ) = {(2,3)} (2.38) 

C 5 := C(e 5 ) = {(1,4)} (2.39) 

For x > and h G L 2 (M 3 ,d 3 fe) we define 

a h (x) = Anx 2 coth(f3x) [ \h(2x, S)| 2 dS, (2.40) 

Js 2 

and for x = we set 

a h (0) = 4vTlima; 2 coth(/3x) / \h(2x, S)| 2 dS. (2.41) 

Furthermore, let 

Y 2 = \Z[4 K 2 4r 2 -li(\l + p 2 ) 2 a 2 g (B 2 ) (2.42) 

y 3 = p^KlKy-l^Xl + plfaliB^-Am^ (X 2 + p 2 ) rr'.f 2 \, (2.43) 

5 Since the density matrix is hermitian, it suffices to know the evolution of the elements on and 
above the diagonal. 
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(principal value square root with branch cut on negative real axis) where 



r = P.V. / j^d 3 k, r;- =4vrW \g(2Bj, £)| 2 d£. (2.44) 

\k\ J S 2 

Theorem 2.3 (Decoherence and thermalization rates) Take coupling functions 
in (2.8) -(2.11) satisfying f 1 = f 2 = f , g\ = g 2 = g. The thermalization and decoher- 
ence rates are given by 



7 th 


= min {(Aj + ^(Bj)} + 0(x 4 ) 


(2.45) 


7 2 dcc 








-y 2 + (^ + Z / 1 2 ) f 7 / (o) + o(x 4 ) 


(2.46) 


7 3 dcC 


= I (A? + fif)o-g(Bi) + \{\ 2 2 + f4)v g (B 2 ) 






-y 3 + (^ + ^ 2 >/(o) + o(x 4 ) 


(2.47) 


7f C 


= (A 2 + ^ 2 ) f r 3 (i? 1 ) + (A 2 + /i 2 K(i?2) 






+ [(ki - k 2 ) 2 + v\ + z/ 2 ] (7/(0) + 0(x 4 ) 


(2.48) 


7 5 dcC 


= (A 2 + ^ 2 )a 3 (i? 1 ) + (A 2 + /x 2 K(S 2 ) 






+ [(Kl + K 2 ) 2 + ^ + v\\ (7/(0) + 0(X 4 ) 


(2.49) 



We give a proof in Section 4. 

Discussion. (1) The thermalization rate depends on energy-exchange couplings 
only. This is natural since energy-conserving interactions leave the populations invari- 
ant. 

(2) For the purely energy-exchange model (/Cj = Uj = 0) the rates depend sym- 
metrically on the local and collective interactions trough the combination A 2 + // 2 . 
However, for the purely energy-conserving interaction (Aj = fj,j = 0) the rates are not 
symmetric in the local and collective interaction. For instance, 74 6C may depend on 
the local interaction only (/ci = k 2 ). 

(3) The effect of energy-conserving and energy-exchange couplings, and of local and 
collective ones, are correlated. Indeed Y 2 , Y3 are complicated functions of all interaction 
parameters, except the local energy-conserving ones (vj). 

2.2.4 Entanglement dynamics 

Consider the family of pure intital states of S given by po = IV'KV'I) with 

V> = 7 |2 \, |2 l++)+ /{ = |2 I--), (2-50) 

where | + +) etc are defined after (2.15), and where 01,02 € C are arbitrary. The 
concurrence is 

C(w)=2 kftkf- (2 - 51> 
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This class of initial states covers the whole range from unentangled states (e.g. a\ = 0) 
to maximally entangled states (e.g. a\ = a 2 £ K). According to Theorem 2.1, the 
density matrix of S at time t > is given by 



Pt 



pi a 

p 2 

p 3 

a P4 



+ 0(x 2 ), 



(2.52) 



with remainder uniform in t, and where pj = Pj(t) and a = a(t) are given by the main 
term on the r.h.s. of (2.24). The initial conditions are 



Pl( Q )= I- 12 4,1 12 ' P2(0) =P3(0) =0, P4 (0) 



|«l| + l a 2| 



|«l| 2 + |fl2| 2 ' 



and 



q(0) 



a\a 2 



|«l I 2 + I cc 2 1 2 ' 



We define 



and note that 



p:=pi(0) G [0,1] 



(2.53) 
(2.54) 

(2.55) 
(2.56) 



p 4 (0) = 1 - p and |a(0)| = ^p(l -p). 

In terms of p, the initial concurrence, (2.51), is C(po) = 2y / p(l — p). There are four 
resonances associated to the eigenvalue e = 0. One is located at the origin, two have 
leading (in x) imaginary parts given by (see (4.10)) 

5 2 := (A 2 + M1KGB1), h := (X 2 2 + p 2 2 )a g (B 2 ), (2.57) 

and the fourth one has leading imaginary part 5 2 +5s. The leading term of the imaginary 
part of e 2 ( Bl+B2) is (see (2.49)) 

65 ■= 5 2 + 5 3 + [(/ei + k 2 ) 2 + v\ + 4] (7/(0). (2.58) 



We also define 



5 + := m&x{5 2 ,5 3 }, 5- := min{<5 2 , £3}. 



(2.59) 
that 



Theorem 2.4 (Disentanglement time bounds.) Take p ^ 0,1 and suppose that 
(52 , <?3 > 0. There is a constant x$ > (independent of p) such that for / x < 
x o VpO- ~~ P) we h ave the following. 

A. (Entanglement death.) There is a constant Ca > (independent ofp,x) s.t 
C{pt) = for all t > tA, where 



1 , 

ia ■= max < — In 

05 



C A 



VpU- -p) 



'^2+^3 



In 



C A 



pO--p) 



C A 

$2 + <^3 



(2.60) 



14 



B. (Entanglement survival.) There is a constant C B > (independent of p, x) 
s.t. C(pt) > for all t <t B , where 

t B := min j^-L- ln[l + C B p(l ~p)},j^ In [l + C B x 2 ] , - ^ ^ } . (2.61) 

Discussion. (1) The disentanglement time is finite since 62,63 > (which implies 
thermalization) . If the system does not thermalize, then entanglement for certain 
initial conditions may stay nonzero for all times. 

(2) The rates 5 are of order x 2 . Both tA and t B increase with decreasing coupling 
strength. This is consistent with the expectation that disentanglent happens at a slower 
pace for small couplings (no disentanglement for x = 0). 

(3) As functions of p, both tA and t B are maximal at p = 1/2 and minimal at 
p £ {0, 1}. This is consistent with the expectation that a maximally entangled state 
(p = 1/2) keeps its entanglement for longest. Even if the initial state is disentangled 
(p = 0, 1) we expect to see creation of entanglement due to the collective coupling 
(up to times at most Ca/{&2 + 63), see (2.60) with p = 0, 1). We show numerically in 
[21] that the resonance dynamics does reveal creation, as well as death and revival of 
entanglement. 

3 Proof of Theorem 2.1 

In [23, 24, 25] we have proven a representation of the reduced dynamics of a general 
iV-level system coupled to a heat reservoir. We outline here how to generalize it to the 
present situation of two spins coupled to three reservoirs. Of course, the generalization 
to a general iV-level system in contact with any number of reservoirs is immediate. 

Resonance energies bifurcate out (x 7^ 0) of the real energies of the system Li- 
ouvillian L§ (2.19), into the upper complex plane, c.f. (2.21). They are the eigenvalues 
of a closed (unbounded, non-normal) operator 

8 

K }it e = L o + 0N + Y l XjW j (3.1) 

i=i 

acting on the GNS Hilbert space 

n Gm = n®'H, (3.2) 

see (2.6). Here, Lq = L$ + + £r 2 + £r 3 with Ls = H$ <8> 1 — 1 <8> Hs and the 
Ljt_. are defined similarly. We write ok for the collection of all the coupling constants 
in (2.8)-(2.11). 6 is a complex spectral translation parameter ([23, 25, 20, 3]) and 
N = Ni + N2 + ^3 is the sum of the total number operators of each reservoir, with 
Nj = N <8>1 + 1<8> N acting on ^J-^, N being the number operator on Fock space 
J 7 . To ease notation, we have labelled in (3.1) the interaction constants by Xj, and 
the Wj are interaction operators, obtained from (2.8)-(2.11) by passing to the GNS 
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space and adding suitable operators in the commutant - making (3.1) an operator of 
the 'C-Liouvillian type'. The defining property of the Wj is that 

K Z fi fis ® ^R!,R 2 ,R = 0, 

where (see also (2.1) and (2.15)) 

! 4 

k=i 

is the trace state of 51 + 52, and where fiR 1; R 2j R ^ s the product state of (six) vacuum 
states fi in each of the three (doubled) factors of the Fock spaces. We refer to [23], 
Appendices A and B for the construction of the explicit expressions of the Wj. The 
following result is the analogue of Theorem 4.1 of [23]. 

Theorem 3.1 (Uncovering of resonances) Suppose that condition (A) is satisfied. 
Fix any 9\ with < 6\ < 9$ (see condition (A)). There is a constant xq (depending 
on 9\) s.t if k := |x| < xq and 9\ < 9$ < 9$, then the operator K%q has only isolated 
eigenvalues in the region {z£C: 9z < #i/2}. These eigenvalues, denoted e{ s \ do not 
depend on 9 and have the expansion 

4 s) = e + # + 0(x 4 ), (3.4) 
where e G spec(Ls), 1 < s < mult(e) and 5^ G C satisfies > 0. Furthermore, 

(s) 

the Se are the eigenvalues of the level shift operator 

A e := -P e WP e (L - e + iO)' 1 P e WP e [ RanPe , (3.5) 

where P e is the spectral projection of Lq associated to the eigenvalue e, P e = 1 — P e , 
and where we have set for short W = X^j=i K j^j- 

Let r be a simple closed contour containing all the eigenvalues of K but not 
containing any of its continuous spectrum, and let 

Q = ^.j^K jifi -z)- 1 &z (3.6) 

be the associated Riesz projection. It has dimension sixteen in our model of Section 
2.2.3, which is the number of eigenvalues inside T, counting multiplicity. The operator 
K$q is reduced by Q and has a finite-dimensional block QK^gQ. We define 

with P R = |fi Rl , R2 , R )(fi Rl ,R 2 , R |, (3.7) 

where the trace is taken over the doubled Fock spaces of all the reservoirs. For each 
t > 0, V(t) is an operator acting on the GNS space Hs^Hs. The trace state fis, (3.3), 
is cyclic and separating for the von Neumann algebra Wis = B(Hs) ®1 C B(Hs ® 'Hs); 
and so the relation 

(V(t)A) ®1Q S = V(t)(A i)fi s VA e B(H S ) (3.8) 

defines uniquely a linear operator V(t) acting on B(Hs) (Heisenberg evolution). 
Theorem 3.1 of [23] can be formulated in the following way. 



V(t) = Tr Rl+R2+R 



Pr e 



itQKsgQ 
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Theorem 3.2 ([23]) Assume the conditions of Theorem 3.1. Then there are constants 
xq, C s.t. if x < xq then 



\ Us (of(A))-us{V(t)A)\ <Cx 2 e-^, 



(3.9) 



for all t > and all A G B{T-L^). Here, 7 = 96* — O(x) satisfies 7 > 29e^ /or a// e, s. 
3.1 Proof of Theorem 2.1 

The dynamical map V(t), (3.7), does not have the group property in general [which 
would mean V(t)V(s) = V(t+s)] and hence nor does V(t). However, due to assumption 
(F) saying that all resonance energies are simple, we have the diagonal representation 



mult(e) 

E E ^ s) \xi s Mxi% 

e s=l 



(3.10) 



where the double sum is over all resonances (see Theorem 3.1), and where K^^x 



(s) 



te^Xe^ and [K^fiYxi' = sf'xe' (adjoint operator), with normalization 

xi s) ,x^) = W 

(Kronecker deltas). An expansion in x yields 

\xi s) )(xh = \vi s) )(vi s) \®PR + d(K), 



(3.11) 



(3.12) 



where the rjf^ and rff^ satisfy (2.20), Pr is defined in (3.7) and where the remainder 
term satisfies PrO(x) = 0(x 2 ). It follows from (3.7), (3.10) and (3.12) that 



mult(e) 



y W=E E e ife ^[|^))(^)|+0(x 2 ) 

e s=l 

Next, we have (recall that the are given in (2.15) and Og in (3.3)) 

(|*n)(*m|®l)fis 

= ^(|^)(^| ® l)fi s ® (#,^ ® # m ) . (3.14) 

k,l 

Combining (3.13) and (3.14) with (3.8), we obtain 
V(t)\$ n ){$ m \ 

mult(e) 

■EE' 



(3.13) 



ile 



e s=l 



E (*l ® **• iP) (^\*n ® * m ) |*,><**| + 0(X 2 



. (3.15) 
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Note that rji^ , rf e s ^ G RanP e (spectral projection of L§ associated to e) and therefore 
the main term of the sum vanishes unless e = E\ — E^ = E n — E m . Thus 

m\) 

va\At{E n -E m ) 

s=l (k,l)£C(E m -E n ) 
XU} S (\$l)($k\) +0(x 2 ) 

A t (m,n;k,l) [ Po ]ki + 0(x 2 ), (3.16) 

(fc,/)eC(£ m -£„) 

where A t is given in (2.27). Combining (3.16) with Theorem 3.2 yields (2.24). 

Finally we check the Chapman-Kolmogorov equation (2.25). Write E; L j for E\ — Ej 
and fyj for <£j ® <Pj. The r.h.s. of (2.25) equals 

mult(E nm ) mu\t(E qp ) , 

X X g &nm J^qp 

(p,q)eC(E mn ) s=l s'=l 

» \ /^) * \ /*.. /«(«') 

/ \ ' E 9P ' 

Since (p, qr) G C(E mn ) we have .E gp = E nm = e, and we can perform the sum over p, q 
in (3.17), 

J2 (4 S '\* qP ) (<?V,#) = (v { /\vi s) ) = S s , s ,, (3.18) 

(p,q)eC(E mn ) 

where we use Yl(n q)eC(E nm ) \ < ^qp)(^qp\ = i n the first step, and (2.20) in the second 
one. Therefore, (3.17) becomes 

mult(£„ m ) 

E e'^£™ (^1,^), (3.19) 

S=l 

which is A t+r (m, n; k, I). Thus (2.25) is proven. It is also easy to establish (2.26). This 
completes the proof of Theorem 2.1. ■ 



x (WeL) (tSL*™) • ( 3 - 17 ) 



4 Level shift operators, proof of Theorem 2.3 
4.1 Level shift operators 

The general form of the level shift operator (defined in (3.5)) with interaction linear in 
creation and annihilation operators has been given in [23], Proposition 5.1. 6 We do not 

6 In the present work we use the trace state for the two spins as reference state (3.3), while in [23] 
the Gibbs state was used - the corresponding modification in the level shift operator is obtained simply 
by setting /3s = 0. Furthermore, the present definition of the level shift operator differs from that of 
[23] by a sign. 
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present the explicit calculations, which are rather standard (albeit a bit lengthy), but 
give the results only. It suffices to consider A e with e > 0, i.e. with e = ej in (2. 34). 7 
For x > and h £ L 2 (M 3 ,d 3 /c) we define 



a±(x) = 2TTX z - 



|/i(2x,£)| 2 d£. 



(4.1) 



sinh(f3x) J S 2 

• The eigenspace of Ls associated with e = is spanned by {<P\ <8) ^2 ® &2, ^3 
^3, ^4 <8> ^4}, where the <Pj are given in (2.15). In this basis we have 

A = iirfa-^ + Xla^B,)} * ' n '[ ' (-1.2) 



+i{l£<T- {B 2 ) + \la-{B 2 )} 
















e 2/3Bi 





_ e 2/3Bi 


-1 





1 








-1 





1 




- e 2/3B 2 


_ e 2/3B 2 










-1 


1 








\ 








e 2/3B 2 


_ e 2^S 2 










-1 


1 



(4.3) 



We see that Ao is the sum of two terms (4.2) and (4.3), representing the (independent) 
interaction of the reservoirs with spin 1 and spin 2, respectively. Ao does not depend on 
the energy-conserving interaction (on xi j2 and vi >2 ). Indeed, those interactions leave 
the diagonal of the density matrix (in the energy basis) invariant. The contributions 
coming from the local (^1,2) and from the collective (Ai^) couplings enter Ao in the 
same way. 

• The eigenspace of Ls associated with e = 2B± is spanned by {^1 <8> #3, #2 <8> ^4}, 
where the <Pj are given in (2.15). In this basis we have 



A2 Sl = {i[/#T-(Bi) + A?<Tj(Bi) 



+ 



filrg^ + XjrgiB,)] 



1 
1 



,2/3^2 



-1 



1 
1 



+i{» 2 2 a-(B 2 ) + \ 2 2 a g -(B 2 )} 

+i[AC 2 a / (0) + ^ 1 V /l (0)] 

— 2k\k 2 t 

where r is given in (2.44) and 

r s (x) = J P.V. / u 2 |<7(|u|,S)| 2 coth 
1 JRxS 2 



_ e 2/3B 2 
1 



1 
-1 



2 



l 



u — 2x 



dudS. 



(4.4) 



(4.5) 



• The eigenspace of Ls associated with e = 2£>2 is spanned by {^1 <g> ^2,^3 <8> ^4}, 
where the <Pj are given in (2.15). We obtain A2B2 in this basis by switching all indices 
1 -f-T- 2 labelling the spins in (4.4). 



7 It is not hard to see that A_ e = — JsA e Js, where Js is the modular conjugation associated to the 
von Neumann algebra B(Us) ®lC B(Hs ® Ws) and the trace state (3.3). 
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• The eigenspace of L$ associated with e = 2(B 2 — B\) is spanned by ^3 <8> ^2> 
where the 0j are given in (2.15). The level shift operator A 2 (£ 2 _£ 1 ) ls J us t a number 
(times the projection operator \<P S <P 2 )(<P S ® ^2!)- We have 



A 



■zau- = + A?cj 9 (Si) + [i\a g2 {B 2 ) + A^cr 9 (S 2 ) 

+ («1 - K2)V/(0) + v\o S M + l^<7/ 2 (0) 

+M?r fll (Bi) + \\r g {B{) - fi 2 2 r g2 (B 2 ) - X 2 2 r g (B 2 ). 



(4.6) 



• The eigenspace of Ls associated with e = 2(B 2 + B\) is spanned by ^1 <8> ^4, 
where the <Pj are given in (2.15). The level shift operator A 2 (£ 2 _|_£ 1 ) is just a number 
(times the projection operator \<Pi <g) ^4)^1 <S> ^4))- We have 



A 



2(B 1 +B 2 ) 



= 1 



lila gi {B 1 ) + \\a g {B l )+ l ila g2 {B 2 ) + \la g {B 2 ) 

+(ki + «2) 2 <7/(0) + ^^(0) + I/f<7fr(0) 
M?r fll (Bi) - Afo(Bi) - M|r fl2 (S 2 ) - A|r fl (B 2 ). 



(4.7) 



4.2 Resonance energies and resonance vectors 



We give here the eigenvalues and eigenvectors of the level shift operators of the last 
section. To ease notation we assume from now on that the form factors governing the 
energy-exchange interactions with all reservoirs are the same, and that those governing 
the energy-conserving ones are too, 



91=92= 9, 



and 



h = h = /■ 



(4.8) 



In order to distinguish different contributions, we still keep all the coupling constans 
distinct. We present here the resonance data: the eigenvalues 5 e and the resonance 
states rji^ and rje\ defined in (2.20). For notational convenience, we set 



3 = 1,2. 



(4.9) 



The resonance vectors below are written in coordinates associated to the bases given 
in the above section. 
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• The resonance data associated to e = are 



c(2) 



0, 



i(A*l + A|)<T ff (S 2 



?( 3 ) - i 



5(4) _ 



§ m , r(3) 
°0 + °0 ' 



(2) 



(3) 



(4) 

Vo 



1 
1 
1 
1 

-e 2 
1 

-e 2 
1 

-ei 
-ei 

1 

1 

eie 2 
-ei 
-e 2 
1 



~(2) 



~(3) 



~(4) 



Tre^S 



Tre-f> H S 



y/fg/f 
Tre"^ 



l/ y/eie 2 
\ Ze 2 /ei 

-1/ei 
1/ei 
-1 
1 

-l/e 2 
-1 

l/e 2 
1 



(4.10) 



1 1 



• The level shift operator (4.4) has the form 
A 2Bl = All + 



e 2 B + C -e 2 B 
-B B-C 



(4.11) 



where A = i(A 2 +/x 2 )±a !? (B 1 ) + i( K 2 + uf)a f (0) - (A? + ^ 2 )r 9 (Bx), B = i(A 2 + // 2 )<r-(B 2 ) 
and C = -2KiK 2 r (recall definitions (2.40), (2.41), (2.44) and (4.1), (4.5), (4.9)). 
The resonance energies associated to e = 2Bi are 



*(±) 

} 2B X 



A + \B{1 + e 2 ) ± ± [B 2 (l + e 2 ) 2 + 4C[B(e 2 - 1) + C]] 



1/2 



(4.12) 



where the square root is the principal branch (with branch cut on the negative real 
axis). The associated resonance eigenvectors are 



„(±) _ 

7 l2B 1 — 



1 

y± 



~(±) _ 1 

^Bi - l+e 2 (y±) 2 



1 

e2V± 



(4.13) 



e 2 B 



where y± is the complex conjugate of y± = 1 + 

Remark. If C = then the eigenvalues (4.12) reduce to A + B(l + e 2 ) and ^4 and 
the resonance vectors have easy expressions too. 

• The resonance data for e = 2B 2 is obtained from that of e = 2Bi by switching 
indices 1h2 labelling spin 1 and spin 2. 

• The resoance energy associated to e = 2(B 2 — B\) is given by (4.6), and we have 

V2(B 2 -B 1 ) = %(B 2 -Bx) = ^3 <8> <2>2- 

• The resoance energy associated to e = 2(B 2 + B\) is given by (4.7), and we have 

V2(B 2 +B 1 ) = »te(B 2 +Bi) = #1 ® #4- 
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4.3 Proof of Theorem 2.3 



In order to obtain the thermalization and decoherence rates, we simply need to calcu- 
late the imaginary parts of the second-order contributions to the resonance energies, 
calculated in the previous section, and invoke Theorem 3.1. Relations (2.45), (2.48) 
and (2.49) are immediate from (4.10), (4.6) and (4.7). Relation (2.46) follows from 
(4.12). This completes the proof of the Theorem. ■ 



5 Proof of Theorem 2.2 

For a density matrix diagonal in the energy basis, p = diag(pi,p2,P3,P4), we have 
D(p) = —2 mm{y/pip4, yJP2Pz} (see 2.30). By the property of return to equilibrium, 
the reduced density matrix of Si + S2 approaches the Gibbs state modulo an error, 

lim p t = p s ^ + 0{x 2 ), (5.1) 

where, in the energy basis, 

PS,, = diagCe-^+^.e-^-^.e-^+^.e-^"*)). 
We have D(ps^) = — 2 Tre i )3gg and so (5.1) implies 

lunD( Pt ) = -2^-^ + 0^). (5.2) 

The concurrence vanishes if D(p t ) < 0. By a Dyson series expansion, one can show that 
the error term in (5.2) is uniform in (3 for (3 < /3q, where /?o < 00 is any fixed number. 
(See also [17] for a more detailed analysis and a better bound.) Therefore, if (3 < /3q, 
then there exists a constant C > (depending only on (3q) s.t. if Tre~^ s < Cx~ 2 then 
the right side of (5.2) is strictly negative. Then the existence of a finite to follows from 
the continuity of t 1 — > D(pt). Next suppose that return to equilibrium takes place at 
rate 7, i.e., that \\p t - p S)/3 + 0(x 2 ) || < Ce"* . Then £(p t ) = £(p S;/3 ) + 0(x 2 ) + 0(e"^) 
and by perturbation theory D(p t ) = —2 ^—^^ +0(x 2 ) + 0(e" 7t ). Standard estimates 
on return to equilibrium show that the remainder term 0(e~ 7 *) is uniform in T varying 
in compacta in (0, 00) [23, 25, 17]. Thus there is a constant d > such that D{p t ) < 
for t > 7~ 1 ln[c / Tre- /3Hs ]. ■ 

6 Proof of Theorem 2.4 

The dynamics of the matrix elements in (2.52) are obtained according to Theorem 2.1. 
For instance, with <£>;.,• = $j ® <J>j ((2.15)) 

pi(t) = pA t (ll;ll) + (l-p)A t (ll;44) 

= Ee ite <' ' {p(*n " *uA S) ) + (**iA S) )} - (6-1) 

s=l 
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where the resonance energies and resonance vectors are given in Section 4.2. Since pt 
is a density matrix, the diagonal elements 

x 3 {t) := \pt\jj, (6.2) 

j = 1,2,3,4, must be non-negative and add up to one. It follows from (2.52) that 

Pj {t) = Xj (t) + 0(x 2 ), j = 1,2,3,4, (6.3) 

and 



Pt 



+ 0(x 2 ). 



(6.4) 



x\ a 
x 2 
x 3 
a X4 

We recall that the remainder term in the previous formula, as well in all that follows, 
is uniform in t > 0. It is sometimes more practical to consider the Xj instead of the 
Pj since the former are known to be non-negative. The spectrum of (2.32), with p 
replaced by pt, (6.4), is 



spec(£ t ) = {x 2 x 3 ,x 2 x 3 , y X ±X4 ± |a|] 2 } + 0(x 2 ). 



(6.5) 



Let D = D(p t ) be the quantity defined in (2.31). In order to calculate D, we need to 
know which of the eigenvalues of £ is the largest one. 

We have the following expressions for x±, . . . ,2:4 (see (6.3) and e.g. (6.1) for x\) 



-I3{B 1 +B 2 ) 



Xi 



X2 



X3 



X4 



1 -e"* 52 )(l - e" M3 ) + 



+p {e-* 52 (e 2 + 1) + e-' 5s (ei + 1) + <T u \e x e 2 - 1)}] + 0{ 



hp{e~ tS2 \ 

e -l3(B 1 -B 2 ) 

+(1 - e~ tSi )e^ 1 {p{e l e 2 - 1) + 1}] + 0(x : 

-j3(-B 1 +B 2 ) 



(6.6) 



(1 _ e ~tS2 )e - i {p(e2 + i) _ i } + (1 _ e-^){-p( ei + 1) + 1} 



(6.7) 



Tre-^s 



(1 _ e -^){- p (e 2 + i) + 1} + (i _ e -^) e ^{ p ( ei + 1) - 1} 



+(1 - e-* <54 )er 1 {p(eie 2 - 1) + 1} + 0(x 2 ) 



e -/3(-B!-B 2 ) 



Tre"^ 
+p(l-e-* <52 )(l-e-*' 53 



(l-p)(l + e 2 - 1 e-* <52 )(l + er 1 e-^) + 
+ 0(x 2 ). 



(6.8) 



(6.9) 



(6.10) 



In the above expressions, we have set for short 

ej ■= e 2 ^ and S s := 5^ , s = 2, 3, 4 

(recall that ^q 1 " 1 = and that 64 = 5 2 + ^3, see after (2.57)). We have also used the 
mean value theorem to obtain the estimate 

e «4 s) = e ^+°(* 4 )) = e 1 ^ + 0(x 2 ), 3 = 2, 3, 4, 
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which holds uniformly in t > 0, provided 5 S > 0. 
The off-diagonal matrix element is estimated by 

a(t) = e ite2 < s i+ s 2>a(0) = e itSR£2 ( s i+ s 2)e-* 55 c»:(0) + 0(x 2 ), (6.11) 

where 5$ is given in (2.58), and where a(0) is linked to 01,02 by (2.54). 

The above expressions for x\, . . . X4 can be used to arrive at the following result. 

Lemma 6.1 (a) We have, uniformly in t > and in p £ [0, 1], 

X1X4 - x 2 x 3 = e^ 62+5:i) p(l-p) + 0(x 2 ). (6.12) 

(b) We have, uniformly in t > 0, 

xix 4 > [Tre _/3 ^ s ]" 2 p(l - p) + 0(x 2 ). (6.13) 

Proof of Lemma 6.1. Relation (6.12) is obtained by direct calculation from (6.6)- 
(6.9). To show (b), we note that by (6.6), and since 1 > p, 

e - / 3(B 1+ B 2 ) 

xi > Tre _^ s pf(t) + 0(^), 
with f(t) = 1 + e' t5 ' 2 e 2 + e - * 53 ei + e'^e^ > 1. Moreover, by (6.9), 

e -0(-Bi-B 2 ) 
X4> Tre .^ g (1- P ) + 0(X 2 ), 

uniformly in t > 0. This ends the proof of Lemma 6.1. ■ 

In what follows we distinguish two cases (inequalities A and B below). We denote 
by C generic constants which are independent of p and x, but whose values may change 
from expression to expression. 

A x 2 x 3 > [y/xjxl + \a\} 2 + 0(x 2 ). 

We first prove that if t > tA (see (2.60)) then inequality A holds. It follows from 
(6.12) that inequality A holds if e -^ 52+5 ^p(l - p) + 2\a\y/x~[xl+ \a\ 2 < Cx 2 , for 
some C > 0. Since Xj < 1 + 0(x 2 ) and \a\ < e~ t5 -' a/p(1 - p) + 0(>c 2 ) (see also 
(6.11) and (2.56)), this condition is satisfied provided 

e^ 5 Vp(l-p) < Cx 2 and e^ (<52+<53 V(l - p) < Cx 2 . (6.14) 

Conditions (6.14) hold if t > tA, see (2.60). Note that only the first two terms in 
the max on the r.h.s. of (2.60) are needed for this argument, the last one will be 
used below. 

Next we show that 

D = -2max{^/xTxI, \a\} +0(x 2 ). (6.15) 
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Inequality A implies that the largest eigenvalue of is x 2 x 3 + 0(x 2 ), see (6.5). 
To calculate D, (2.31), we need to take the square roots of the eigenvalues of £. 
Using (6.13) together with (6.12) we obtain 

x 2 x 3 > Cp(l -p)- e- t{52+5s) p(l -p)+ 0(x 2 ), (6.16) 

and consequently, since t > t a and hence t > C s ^ s , we have 

x 2 x 3 >Cp(l-p) + 0(x 2 ). (6.17) 



We conclude that for p ^_ p ^ small enough (|x| < XQy/p(l — p) for some xq in- 
dependent of p) , we have the following expressions for the square roots of the 
eigenvalues of 



\y/x^xl± \a\) 2 + 0{x 2 ) = \y/x^xl± \a\\ + 0(x 2 ) (6.18) 
y/x 2 x 3 + 0{x 2 ) = ^ii + O^ 2 ). (6.19) 

Using expressions (6.18) and (6.19) in (6.5) we arrive at (6.15). 

We are now ready to complete the proof of point A of Theorem 2.4. We have 
\a\ < Cx 2 , see (6.14), and hence 

max{^/xiX4, \a\} > Cy/p(l — p). 

Therefore, by (6.15), D < and by (2.30), the concurrence vanishes. 

B x 2 x 3 < [y/xixi + \a\} 2 + 0{x 2 ). 

Due to (6.12) and (6.13), inequality B holds if 

e- t{52+53) p(l-p) + Ce~ t5r 'p{l-p) + e -2 *S(l - p) > Cx 2 , 

for some C > 0. The latter condition is satisfied if either of the tree positive 
summands on the left hand side are bounded below by Cx 2 , i.e. if 

t < m a JUn(c B P{1 - p) ),^ln(c B P{1 - p) 



$5 V ) ' 62 + 63 V k2 

1 Infc^-^V (6.20) 



§1 + h \ 

We have used here that > 5 2 + S3, see (2.58). 

Next we analyze D, (2.31). The largest eigenvalue of £t, (6.5), is v\ = \yjx\x± + 
|a|] 2 + 0(x 2 ). Its square root is y/u 1 = y/xix^ + |a| + 0(x 2 ) (this follows from 
(6.13) and \x\ < x§\Jp(\ — p)). The quantity D is obtained by subtracting from 
the terms yj (yjx\x~l — \a\) 2 + 0{x 2 ) and twice y x 2 x 3 + 0{x 2 ). We are now 
showing that for t < t B , we have 

Vx 2 x 3 + 0(x 2 ) < ^v^i and (6.21) 

^(^l-\a\) 2 + 0(x 2 ) < 1^. (6.22) 
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It then follows that for t < ts, we have C(pt) > joy/vi, an d due to (6.13) the 
statement of point B in Theorem 2.4 holds. 

It remains to show (6.21) and (6.22). We obtain the upper bound X2X3 < C(l — 
e _t<54 ) 2 +0(x 2 ) directly from expressions (6.7), (6.8). By taking this into account, 
together with (6.12), (6.13) and |x| < xoy / p(l — p), we see that (the square of) 
(6.21) holds provided that 1 — e~ t6i < Cp{l — p) (where C is small), which in 
turn is implied by 

t < — \n[Cp{l - p) + 1]. (6.23) 
64 

To summarize, condition (6.23) implies (6.21). 

Our next task is to prove (6.22). By squaring this inequality we see that it is 
satisfied provided that (^x\X4 — \a\) 2 < 1^X1X4 + 0{x 2 ). Invoking (6.12) and 
|x| < xq^/p{1 — p), we see that the last inequality holds provided 

(v 7 ^- M) 2 < Cp{\-p). (6.24) 

Note that the l.h.s vanishes at t = 0, so the inequality holds for small enough 
times. Let us set 

5- := min{c)2, 63} and 5 + := max{e>2, ^3}- (6.25) 

It follows from (6.6) and (6.9) that 

X1X4 < [TYe-^ s ]- 2 |p(l -p)e- ts -( ei + e 2 + e x e 2 + 1) 

x(l + e- t<52 / ei )(l + e-^/d) + C(l - e-* 5 +) 2 } + 0(x 2 ). (6.26) 

We estimate 

(1 + e- t<52 / ei )(l + e- t<5 Vei) < (1 + l/ei)(l + l/e 2 ) + C(l - e"* 5 +). (6.27) 

Combining (6.26) and (6.27), and using the definition (6.10) of ej, we obtain the 
upper bound 

X1X4 < p(l - p)e~ tS - + C(l - e' t5 +) + 0(x 2 ). (6.28) 

Furthermore, if t < ^j- ln[l / (1 — Cm 2 )} for some constant C > 0, then C(l — 
e~ tS +) = 0{k 2 ). The last upper bound on t is implied by 

t < ^-\n[l + Cx 2 }. (6.29) 

Next we have |a| = e~ tSs \fp{l — p) + 0(x 2 ), so we obtain for t satisfying (6.29) 
V^ixl - \a\ < y/p(l - p) \e'^ s - - e ~ t55 ] + 0(x 2 ). (6.30) 
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We also know that — Cx 2 < yjx\x& — \a\, a fact that follows simply from the 
positivity of pt (see also (6.4)). Therefore, 



(s/xiXi — \a\) 2 < max < Cx 4 , p(l — p) 



+ 0(x 2 ) 



We use this upper bound to see that (6.24) holds provided 

1 2 



p(l-p) 



,-t<5 5 



< Cp{l -p) + 0(x 2 



(6.31) 



(6.32) 



Since |x| < x,§\Jp(\ — p) inequality (6.32) is implied by the bound e 



-tS B 



< 



C. We have e - ^- - e~ tSs = e - ^- 
holds if 1 - e -t(S 6 -S-/2) < Cf which 

in turn is implied by 

C 



[1 _ e -t(* 5 -«-/2)] < 1 _ e -t(« B -«_/2) ) so (6 31) 



t < 



5 5 -S-/2' 



(6.33) 



We have thus shown that if t satisfies (6.20), (6.23), (6.29) and (6.33), then the 
bounds (6.21) and (6.22) hold. Condition (6.20) is implied by condition (6.23) 
since |x| < xoy / p(l — p). Thus all above constraints on t are verified for t < ts, 
see (2.61). This shows point B of Theorem 2.4. 



The proof of Theorem 2.4 is complete. 



7 Davies generator and level shift operators 

We take all the coupling constants in (2.8)-(2.11) to be proportional to some x. Let 
a4 denote the reduced Schrodinger dynamics of the qubits. 

Proposition 7.1 Let p be any density matrix of the qubits. Then 

lim sup |K T/ " 2 o </- 2 (p) - e TK *p\\ = 0, (7.1) 

*-»-0 r>0 

where the operator K* : Z 1 (7^s) ~~ ^ ^C^s) maps density matrices into density matrices. 
Viewing I 1 (H s ) asHs®H s , is an operator on the latter tensor product. As such, it 
leaves spectral subspaces of L%, (2.19) invariant. On the subspace RanPi s=e , if" acts 
as (iA e )* (adjoint of level shift operator (3.5)). 

Remarks. (1) Relation (7.1) implies that e rK> maps density matrices to density 

matrices: indeed, a oaj (p) is a density matrix. if" is thus the Davies generator. 

(2) As defined in (3.5), A e is proportional to x 2 , but in the proposition, we prefer 
to reinterpret the level shift operator to be independent of x (divide (3.5) by x 2 ). 

Proof of Proposition (7.1). Our proof is partly inspired by [13]. We have 

W ° <{p)\mn = <^ ,e i * i -e- i * L «(|1>n)(1> m | ® Ir)^^^"^^) 

= e -itE mn Mm,n;k,l)[p] kl + 0(>c 2 ). (7.2) 

■n/rn, } 
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Here, V&o is the initial state represented in the GNS space and is the standard 
Liouville operator implementing the dynamics. We use that Lq = Lg + Lr, Lg = 
H s (g>l- l(g> H s , H s ® n = E n $ n , and Theorem 2.1. We denote E mn = E m — E n and, 
below, <& m n = &m <8> *n- Setting i = r/x 2 , using the explicit form of A t , (2.27), and 
the fact that the remainder in (7.2) is uniform in t, we obtain (7.1) with 

[e TKi p]mn = Yl Al(m,n;k,l)[p] kl (7.3) 
(k,i)eC(E mn ) 

4(m,n;fc,Z) = (* mn , e - iT ( A ^»)'$ w ) . (7.4) 
To arrive at (7.4), we make use of the diagonalization formula 

mult(B nm ) 

Z / 1 '&nm' N 'ftnm 1 

S=l 

implying (see (7.2) and (2.27)) 

Al(m,n;k,l) = <^,e irA ^$ nm ) 

= J s e irA ^ m J s $ mn ) 

= ($mn,Jse- iT{AE ™ r J S $kl), (7.5) 

where Js is the modular conjugation associated to the pair (9Jts, ^s)- Finally, one sees 
readily (see also [22]) that JsA e Js = — A_ e , so that (7.4) follows from (7.5). 

As an operator on ^(Hs), has matrix elements K mn ^i defined by K$\&k)(&i\ = 
T, m ,n K mn,ki\$ni)(<S , n\, and by applying d T \ T=0 to (7.3), (7.4), one obtains 

This completes the proof of Proposition (7.1). ■ 
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